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We analyse the clustering of cosmic large scale structure using a consistent modified gravity 
perturbation theory, accounting for anisotropic effects along and transverse to the line of sight. The 
growth factor has a particular scale dependence in f(R) gravity and we fit for the shape parameter 
fno simultaneously with the distance and the large scale (general relativity) limit of the growth 
function. Using more than 690,000 galaxies in the Baryon Oscillation Spectroscopy Survey Data 
Release 11, we find no evidence for extra scale dependence, with the 95% confidence upper limit 
|/ho| < 8 x 1(T 4 . Future clustering data, such as from the Dark Energy Spectroscopic Instrument, 
can use this consistent methodology to impose tighter constraints. 

PACS numbers: 98.80.-k;04.50.Kd;98.65.Dx 


I. INTRODUCTION 

The growth of large scale structure in the universe is 
a multifaceted probe of cosmology. The galaxy cluster¬ 
ing pattern measures cosmic geometry through baryon 
acoustic oscillations, giving angular distances transverse 
to the line of sight and radial distance intervals, or the 
Hubble parameter, along the line of sight, and the ratio 
of the two known as the Alcock-Paczyhski effect [1]. The 
evolution of the clustering amplitude provides the growth 
factor and the shape of the clustering correlation function 
or power spectrum depends on early universe conditions 
and later, scale dependent effects. In addition, redshift 
space distortions (RSD) cause anisotropy in the cluster¬ 
ing between the transverse and radial directions, and this 
probes the velocity field and the law of gravity [2-7]. 

General relativity (GR) predicts scale independent 
growth in the linear perturbation regime and specific 
redshift distortion patterns, and so probing for scale de¬ 
pendence or distortion deviations can test the theory of 
gravity. Our aim is to investigate the general relativ¬ 
ity cosmological framework by fitting the galaxy cluster¬ 
ing data while allowing for scale dependence, and con¬ 
strain such deviations. We are particularly motivated by 
scalar-tensor theories and use a perturbation theory tern- 
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plate derived for f(R) gravity. Scale dependence arises 
at length scales smaller than or of order the inverse of 
the scalaron mass m = where a subscript R 

denotes a derivative with respect to the Ricci scalar [8]. 
For GR, / = R and so the scale dependence vanishes. 
This scale dependence was tested using the combination 
of cosmic microwave background (CMB) data and galaxy 
clustering spectra in the following work [9]. 

By employing the perturbation theory of Taruya et 
al. [10] that uses a resummed propagator to partially in¬ 
clude nonperturbative and screening effects of the mod¬ 
ified gravity, we can analyze the clustering correlation 
function to smaller scales than linear theory or simple 
perturbation calculations. We join this to our previ¬ 
ous, substantially model independent approach of treat¬ 
ing the background expansion in terms of the angular 
diameter distance and Hubble parameter [11-13], rather 
than assuming a dark energy model such as ACDM. Fur¬ 
thermore, we generalize the previous scale independent 
growth factor to two quantities: one scale independent 
(corresponding to the large scale limit of the growth fac¬ 
tor) and one scale dependent (which can be thought of 
as characterizing the scalar-tensor modification, e.g. m 
or /ho)- 

While we concentrate here on improving the RSD and 
galaxy clustering analysis for testing GR, various other 
methods have also been used. For example, the Planck 
collaboration has put constraints on a wide of range of 
modified gravity models using the Planck 2015 CMB 
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data, combined with large scale structure observations 
[14]; the SDSS-III (BOSS) team has tested GR using the 
observed structure growth patterns [15-17]; the Wiggle-z 
team has tested modified gravity models [18, 19]; and the 
CFHTLenS team employs the complementarity between 
weak gravitational lensing and RSD [20] . The abundance 
of clusters [21, 22] and the cluster profiles [23, 24] have 
also been used for gravity tests. For more recent obser¬ 
vational tests of GR, see [25-32] and [33] for a review. 

To focus on exploring the scale dependence, we use 
only the clustering data, from the Baryon Oscillation 
Spectroscopic Survey (BOSS) Data Release 11 (DR11) 
of the Sloan Digital Sky Survey 3 (SDSS3) [34, 35]. 
This consists of about 690,000 galaxies over an effective 
volume of 6 Gpc 3 with an effective redshift z = 0.57. 
We measure the two dimensional anisotropic correlation 
function as a function of transverse and radial separation 
between galaxies, and fit this to the redshift space dis¬ 
torted resummed perturbation theory. This comparison 
then imposes constraints on the cosmological and gravi¬ 
tational quantities, allowing us to test general relativity. 

Section II describes the modified gravity theoretical 
approach, from the f(R) gravity model to the resummed 
propagator and resulting perturbation theory to the pre¬ 
diction for the clustering correlation function. We also 
lay out our approach of splitting the growth function into 
scale independent and dependent parts. In Sec. Ill A we 
discuss measurement of the anisotropic correlation func¬ 
tion from the data and treatment of the covariance ma¬ 
trix. We verify in Sec. Ill B that we recover ACDM from 


ACDM simulated data. The comparison of the theory 
to the measurement is in Sec. IIIC and HID, where we 
analyze the fits to the cosmological quantities and their 
consistency with general relativity and ACDM. We sum¬ 
marize and conclude in Sec. IV. 


II. THEORETICAL MODEL 

A. ,f(R) gravity model 

We consider perturbations around the Friedman- 
Robertson-Walker universe described by the metric 

ds 2 = —(1 + 2 T)c!t 2 + a(t) 2 ( 1 — 24>)<5y dx l dx :j . (1) 

We will investigate modified gravity models that can be 
modelled by Brans-Dicke gravity on subhorizon scales. 
The metric perturbations and the scalar field perturba¬ 
tion <f> = (f>o + ip obey the following equations in Fourier 
space [36] 

—/c 2 T = AirGa 2 p5 + — k 2 ip, (2) 

(3 + 2uj B n)\k 2 (p = 8nGpm8 (3) 

a z 

$ - 4/ = tp, (4) 

where X represents the self-interaction, which can be ex¬ 
panded as [36] 


1 f d 3 l?i d 3 k o 

Z(<p) = M 1 (k)ip(k) + - J -^-^—5 D (k - k 1 2)M 2 (k 1 ,k2)ip(k 1 )(p(k 2 ) 

+ ^ J d fc3 Mfc ~ ki 23 )M 3 (k 1 ,k2,k3)tp(k 1 )(p(k 2 )<p(k 3 ), (5) 


where kij = ki+kj and kijk = ki+kj+kk■ We treat the 
matter fluctuations S as a pressureless fluid flow, whose 
evolution equation is given by 

^ + ±V.[(l + 5)v] = 0, (6) 

at a 

<9v 1 1 

—+ J?v+-(v-V)-v =—V*. (7) 

at a a 

We will assume the irrotationality of fluid quantities and 
express the velocity field in terms of the velocity diver¬ 
gence 0 = V • v/{aH). 

In order to develop the template for the redshift space 
power spectrum, we need to specify the interaction term 
T. We consider f(R) models as a representative class of 
models where the linear growth function is scale depen¬ 


dent. f(R) gravity models are described by the action 

S = J d A Xy/^gf(R) + J d 4 x^gC m - ( 8 ) 

We consider the function /(i?) given by the lowest order 
expansion in the small quantity \fm\ 1, 

f{R) = —2k 2 pa + |/flo|-p, (9) 

general to many f(R) models that are observationally 
viable. Here p\ is the constant energy density and Rq 
is the background curvature at present time. For small 
|/ro|j the background expansion in this Ansatz can be 
approximated as the one in the ACDM model. 

The scalar field perturbation is given by 

= Ir — [r 


( 10 ) 
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FIG. 1: ( Left panel) The linear power spectra Pss(k), Pse(k) and Pe©(fc) from top to bottom. The power spectra of the 
ACDM model are presented as black solid curves, and the power spectra of f(R) gravity models with \/ro\ = 3.2 x 10 -5 
and |/ro| = 3.0 x 10 -4 are presented as blue dashed and dotted curves, respectively. The results are evaluated at redshift 
2 = 0.57. ( Right panel ) The growth functions D s + and D ® have a scale independent amplitude, depending on 12 a (taking fixed 
primordial amplitude), and scale dependent shape, depending on \fm\- The results in each subpanel are shown at z = 0.57 for 
|/.Ro| = 3.2 x 10 —5 , and 12 a = 0.68 (black solid curves) and 0.71 (blue dashed curves). 


where /r = df /dR and the bar indicates that the quan¬ 
tity is evaluated on the background. f(R) gravity models 
are equivalent to Brans-Dicke gravity with wbd = 0 and 
the interaction term X is given by 


X = SR=R(f R )-R(f R ). 


Thus the coupling functions M n are 


M n 


d n R(f R ) 

df n R 


( 11 ) 


( 12 ) 


so Eq. (9) gives 


Combining this with the energy momentum conserva¬ 
tion equations, we obtain the equation that determines 
the linear growth factor D+(k,t), 


CD 5 + = 0 , 


C dt 2+2H dt 4nG 


(4 + M ia 2 /k 2 \ 

\z + M 1 a 2 /k 2 ) Pm ' 

(15) 


The growth function D s + and the growth rate D® = 
dD & + /d log a are also scale dependent. We introduce the 
following parametrisation of the growth rates 


D* + (k,t) = G s (t)F s (k,t-,M 1 ), 

D+{k, t) = GeWfefttiMr), (16) 


Mi 


11 R 3 

' 


(13) 


By linearising the evolution equations, the solution for 
the gravitation potential is given by 


fc 2 !? = — 47rG 


/ 4 + Mia 2 /k 2 \ 
\3 + Mia 2 /k 2 ) 


Pmd , 


(14) 


On scales larger than the Compton wavelength of the 
scalar field, m _1 = (3/Mi) 1 / 2 = V3 /rr, the scalar 
field does not propagate and we recover ACDM. On 
the other hand, on small scales, gravity is enhanced: 
G eff (4/3)G. 


where we defined Fg(k,t; M\) and F©(fc,t;Mi) so that 
D s + (k,t) —> Gg(t ) and D+(k,t) Ge(t ) in the limit of 
k —> 0. Since we recover ACDM in the k —> 0 limit, 
Gg and Gq are determined by the usual cosmological 
parameters and they are independent of \/ro\- On the 
other hand, the scale dependence is controlled by Mi in 
Eq. (14), which is determined by \/ro\ and the cosmolog¬ 
ical parameters. 

We will find solutions for S and 9 by solving Eqs. (2), 
(3), (6), (7) using perturbation theory. Once the nonlin¬ 
earity becomes important, the non-linear self-interactions 
M,>i will suppress the scalar field interactions by the 
chameleon mechanism. This effect is included perturba- 
tively in our approach. 
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B. RSD model for f(R) gravity models 

The anisotropy of galaxy clustering is now recognized 
as a useful probe of gravity on cosmological scales. This 
is because the anisotropic clustering signal contains in¬ 
formation on both the cosmic expansion and growth of 
structure, through the Alcock-Paczynski effect and RSD. 
In principle, these two effects are simply described by the 
mapping formula from the statistically isotropic frame, 
however modeling the RSD effect is rather complex be¬ 
cause of the nonlinear and stochastic nature of the map¬ 
ping. As a result, the applicable range of the linear the¬ 
ory prediction is quite limited. Even at the largest scales 
accessible by future galaxy surveys, a proper account of 
the nonlinearity is crucial for a robust test of gravity. 

Here, we will adopt an improved model of RSD by 
Ref. [37] for the theoretical template of the redshift-space 
correlation function. While this model has been orig¬ 
inally proposed to characterize the matter power spec¬ 
trum in GR, the assumptions and propositions behind 
the model prescription do not rely on any specific gravi¬ 
tational theory. Thus it can apply to any model of mod¬ 
ified gravity. Indeed, the model has been tested against 
the dark matter simulation of the f(R) gravity model, 
where a good agreement with AI-body results was found 
[38]. One important remark is that the dynamics of den¬ 
sity and velocity fields in modified gravity can be differ¬ 
ent from GR, and each building block in the RSD model 
needs to be carefully computed, taking a proper account 
of the modification of (non)linear gravitational growth, 
which we will describe below (see also Appendix A). This 
is what we refer to as a consistent analysis of modified 
gravity. 

Employing the linear bias prescription, the improved 
model of RSD is given in Fourier space as function of 
wavenumber k and directional cosine p, = k z /k with k z 
being line-of-sight component of k: 

P{k, p) = {b 2 Pgg{k) + 2[i 2 bGePse(k) + p^G^Peeik) 
+A(k, /x; b , Ge) + B(k, /i; b , G©)} 
xD FoG {kpa p ), (17) 

where b is the linear bias parameter. (See Sec. Ill C for 
a more sophisticated treatment.) The G© is the param¬ 
eter characterizing the growth of structure introduced in 
Eq. (16). The functions Pss, Pee an d Pge are, respec¬ 
tively, the auto-power spectra of density and velocity- 
divergence fields, and their cross-power spectrum. The 
function D FoG characterizes the suppression of the power 
spectrum due to the virialized random motion of galax¬ 
ies [39-41], for which we assume the Gaussian form: 

Dfog{x) = e~ x . (18) 

Since the suppression of the power spectrum basically 
comes from the galaxies sitting in a halo, we shall treat 
a p in the damping function D FoG as a free parameter. 

In Eq. (17), the main characteristic is the A(k , p) and 
B(k,p) terms, which represent the higher-order coupling 


between density and velocity fields. These have been 
derived on the basis of the low-fc expansion from the ex¬ 
act expression for the redshift-space power spectrum, ex¬ 
pressed as: 


3 2 / \ i 

^ = K s EE^”(% 


n— 1 a, 6=1 


k 3 f oo ,1 

dr / dx 


B = b 


(2tt) 2 J 0 j_ x 
x { A ab( r ; x) B 2 ab{p , k-p,-k) 

+Aa b {r,x)B 2 ab{k-p 1 p,-k)Y (19) 


n— 1 a,6=1 


a +k 3 roc n 1 

0 ! Wl, * 1 dx 


'-1 


„ , v Pa2(kVl + r*-2rx)P b2 (kr) 
xB ab {r,x) (1 + r 2 _ 2ra; )a • l 20 ) 


where r = p/k and x = k ■ p/(kp). Here, the func¬ 
tions P ab and B abc are the power spectrum and bispec¬ 
trum of the two-component multiplet ’JG = (<5, 0). The 
non-vanishing coefficients, A" h , A” b , and B™ b , are those 
presented in Sec. III-B2 of Ref. [42] and Appendix B of 
Ref. [37], respectively. 

In calculation of the redshift-space power spectrum, we 
need to properly take into account the effect of nonlinear 
gravitational evolution in each term of Eq. (17). Since 
the standard perturbation theory (PT) [43] is known to 
produce an ill-behaved expansion leading to unwanted 
UV behaviour, we shall apply the resummed perturba¬ 
tion theory called RegPT, which has been formulated in 
Ref. [44] and been later extended in Ref. [10] to the modi¬ 
fied gravity models. Following the prescription described 
in [10], we compute the power spectra Pxv(k) as well as 
the A and B terms, including consistently the nonlinear 
corrections up to the one-loop order in the f(R) grav¬ 
ity model. The explicit expressions for statistical quan¬ 
tities, necessary for a consistent one-loop calculation of 
the redshift-space power spectrum, are summarized in 
Appendix A. 

In the cosmological analysis in Sec. Ill, we will con¬ 
sider the following two cases, which can be thought of as 
fixing the expansion and large scale growth histories, or 
allowing them to vary. In the first case, we will fix all 
the cosmological parameters in the template by a fidu¬ 
cial cosmology and vary \fno\, the linear bias b and the 
velocity dispersion a p . This corresponds to fixing Gg, 
Ge and the amplitude of the initial spectrum to those 
determined by the fiducial cosmology. 

In the second case, we allow G© to float by varying 
Ha = Pa/3H 2 where Hq is the present day Hubble con¬ 
stant in the flat universe. Note that with the variation 
of Ha, the amplitude of the growth function D s + (Gg) as 
well as the shape of the growth function (Fg) and growth 
rate ( Fq ) are also changed through the mass parame¬ 
ter Mi [see Eq. (15)]. In practice, however, we see from 
Fig. 1 that the dependence of the shapes Fg and Fq on 
Ha are very weak. Further, the change of Gg is degener¬ 
ate with the linear bias b in the linear regime. Although 
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FIG. 2: (Left panel) The correlation function £(<7, n) of the ACDM model is shown as black solid contours. Blue dashed 
and dotted contours represent £(< 7 , n) of the ACDM model in which D s + is replaced by that of f(R) gravity models with 
|/.ro| = 3.2 x 10 -5 and \fm\ = 3.0 x IIP 4 , respectively. The thick solid circle represents the BAO ring. The levels of contours 
are given by (—0.005, —0.0025, —0.004, 0.003, 0.005, 0.016, 0.05) from outer to inner contours. (Right panel) Black solid contours 
represent £(cr, 7r) of the ACDM model, and blue dashed and dotted contours represent £(< 7 ,7r) of the ACDM model in which 
Df is replaced by that of f(R) gravity models of \/ro\ = 3.2 x 10~ 5 and |/i?o| = 3.0 x 10 —4 , respectively. 


the one-loop terms break the degeneracy, this effect is 
small in the regime of our interest. Hence, in the second 
case, we vary the scale independent growth rate G@ and 
the linear bias in the PT template, and introduce a new 
parameter 

G h {t) = bG s {t), (21) 

to represent the combined effect of the variation of b and 
G$. Then G$ and the initial amplitude of the power 
spectrum are formally held fixed. 


C. Correlation function £(<7,77) in f(R) gravity 

The two-point correlation function of galaxy clustering 
in the redshift space, £, is described by a function of a 
and 7 r, where a and 7r are the transverse and the radial 
directions with respect to the observer. From the power 
spectrum P(fc,/x), we can compute the correlation func¬ 
tion £(cr, 7 t) by Fourier transformation. The correlation 
function is generally expanded as 

/ r1 3 k 

= E ( 22 ) 

^:even 


with Ve. being the Legendre polynomials. Here, we de¬ 
fined v = 7 t/s and s = (a 2 + 7 r 2 ) 1 / 2 . The moments of 
the correlation function, £f(s), are defined in [11], The 
contributions from moments higher than £ = 8 will be 
ignored in our analysis because we are only interested in 
quasi-nonlincar scales where perturbation theory is ap¬ 
plicable. 

We compute theoretical templates for the moments 
£<(s) using the power spectrum given by Eq. (17). The 
shape of the initial spectra and the growth functions, 
Gs and Ge, are given by the best fit ACDM model 
from Planck 2013 [45]. We provide multiple templates 
with various \fm\ = (0,10 —6 ,3.2 x 10 -6 ,10 -5 ,3.2 x 
10" 5 ,10~ 4 , 3.2 x 10" 4 ,10" 3 ,3.2 x 10" 3 ). For other values 
°f | fno |, we interpolate these templates. 

We now study the variation of £(cr, ir) due to the change 
of the growth function and growth rate. This variation 
was studied in the case of scale independent growth func¬ 
tions in [11] by varying G$(t) and Ge(t) in the RegPT 
template for ACDM. Following this approach, we first 
study the impact of having a scale dependent growth 
function or growth rate. For this purpose, we compute 
the correlation function £ (cr, 7r) using a ACDM template 
and replace the growth function D s + or growth rate D® 
by that in f(R) gravity with |/^o| = 3.2 x 10 -5 and 
\f m \ = 3.0 x 10” 4 . 

For the scale dependent growth function D i, the vari- 
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ation of £(< 7 , n) with a small \f R o\ = 3.2 x 10~ 5 is similar 
to the case of a scale independent enhancement of the 
growth function studied in [11]. Peak points on the BAO 
ring represented by a thick black solid curve in Fig. 2 
move coherently along the circle in an anti-clockwise di¬ 
rection. The blue dashed contours in the left panel of 
Fig. 2 represent this variation. However, £(<r, tt) with a 
larger \f R o\ = 3.0 x 10 -4 varies differently from the scale 
independent case. Peak points on the BAO ring remain 
the same, while minima of BAO are deepened, shown as 
blue dotted contours in the same panel. 

Next, we consider the variation of £(cr, n) due to the 
scale dependent growth rate £>?. In the case of the 
scale independent growth rate, if Gq increases or de¬ 
creases, the anisotropic effects from higher order mo¬ 
ments are visible in the plot of £(<r, n). The location 
of the crossing points between the contour levels and the 
BAO ring (thick solid) shifts clockwise or anti-clockwise 
slightly. The blue dashed contours in the right panel 
of Fig. 2 represent the variation of £(cr, n) with AD® 
for \f R0 \ = 3.2 x 10- 5 and \f R0 \ = 3.0 x 10~ 4 . For 
| f R o | = 3.0 x 10“ 4 , we can see that the peak positions 
are ‘squeezed’ along the BAO ring. 

Having shown the individual effects of a scale depen¬ 
dent growth function and growth rate on the correlation 
function, we now present the correlation function £(cr, n) 
in f(R) gravity models. In Fig. 3, the correlation function 
with |//jo| = 3.2 x 10 -5 and \fm\ = 3.0 x 10~ 4 are plotted 
as black dashed and black dotted contours, respectively. 
There is no variation of £(cr, 7r) up to \f R o\ < 10~ 6 , and 
the correlation function is effectively equivalent to that 
of ACDM. When |/jj 0 | increases to \f R o\ ~ 10 -4 , we 
observe the deviation of £ (cr, n) from ACDM and this de¬ 
viation can be understood as the combined effect of the 
scale dependent growth function and growth rate shown 
in Fig. 2. 


III. METHODOLOGY AND RESULTS 

The observed clustering of galaxies in redshift space 
not only probes the density and velocity fields, i.e. the 
growth and gravity as discussed in the previous section, 
but also provides a useful tool to determine both the 
transverse and radial distances by exploiting the Alcock- 
Paczynski effect and the BAO scale. In galaxy redshift 
surveys, each galaxy is located by its angular coordinates 
and redshift. However, the correlation function, £(<r, 7r), 
is measured in comoving distances. Therefore a fiducial 
cosmological model is required for conversion into comov¬ 
ing space. We use the best fit ACDM universe to Planck 
2013 data. The conversion depends on the transverse 
and radial distances involving Da and H~ 1 . Instead of 
recreating the measured correlation function in comoving 
distances for each different model, we create the fiducial 
maps from the theoretical correlation function by rescal¬ 
ing the transverse and radial distances using Da and H~ 1 
and fit them to the observed correlation function. There- 



50 100 150 200 
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FIG. 3: The best fit correlation function £(cr, 7r) of ACDM 
(black solid unfilled contours) and the correlation function of 
f(R) gravity models with \fm\ = 3.2 x 1CD 5 (black dashed 
unfilled contours) and 3.0 x 10 -4 (dotted unfilled contours). 
The blue filled contours represent the measured £(cr, n) from 
the DR11 CMASS data. The levels of contours are given 
by (—0.001, 0.002, 0.005, 0.016, 0.05) from the outer to inner 
contours. 

fore, when we fit the measured £(er, 7r), the two distance 
parameters of (Da, H~ x ) are added to the structure for¬ 
mation parameter set of {Ga, Go, a p , \f R o\, cr p } discussed 
in Sec. IIB. 


A. Measured £(<r, 7r) using DR11 

Our measurements are based on those previously pre¬ 
sented in [13] which follows a similar procedure to [12]. 

Briefly, in our analysis we utilise data release DR11 of 
the Baryon Oscillation Spectroscopic Survey [BOSS; 46- 
48] which is part of the larger Sloan Digital Sky Survey 
[SDSS; 49, 50] program. From DR11 we focus our anal¬ 
ysis on the Constant Stellar Mass Sample (CMASS) [51], 
which contains 690,826 galaxies and covers the redshift 
range 2 : = 0.43 — 0.7 over a sky area of ^8,500 square 
degrees with an effective volume of V e g ~ 6.0 Gpc 3 . The 
CMASS galaxy sample is composed primarily of bright, 
central galaxies, resulting in a highly biased (b ~ 2) se¬ 
lection of mass tracers [52], 

The redshift-space two-dimensional correlation func¬ 
tion £(cr, 7r) of the BOSS DR11 galaxies was computed 
using the standard Landy-Szalay estimator [53]. In the 
computation of this estimator we used a random point 
catalogue that constitutes an unclustered but observa- 
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FIG. 4: The measured constraints on fno , and their robustness to various tests, are presented. The measured likelihood function 
appears in the top panels and the measured difference of x 2 is i n the bottom panels. (Left panel) Results marginalizing over 
the scale independent growth rate Ge are shown by the black solid curve, while the constraints fixing Ge = 0.46, given by 
the Planck concordance ACDM model, are blue dashed curves. The results for fno do not depend appreciably on the scale 
independent behavior. (Right panel) The results also do not depend significantly on whether the initial power spectrum P(k) 
used matches the Planck (black solid) or WMAP9 (black dashed) model. The blue dotted curve represents the results from 
analyzing galaxy clustering from ACDM mock catalogues, verifying that \fm\ —> 0 is recovered in this case. 


tionally representative sample of the BOSS CMASS sur¬ 
vey and contains ~ 50 times as many randoms as we have 
galaxies. 

The covariance matrix was obtained from 600 mock 
catalogues based on second-order Lagrangian perturba¬ 
tion theory (2LPT) [54, 55]. The mocks reproduce 
the same survey geometry and number density as the 
CMASS galaxy sample. We obtain the covariance ma¬ 
trix using the same treatment presented in our previous 
works [12, 13]. 

We calculate the correlation function in 225 bins 
spaced by 10 h -1 Mpc in the range 0 < ct, tt < 

150 h -1 Mpc. However, at small scales, if the non- 
perturbative effect of FoG is underestimated, then the 
residual squeezing can be misinterpreted as a variation 
in Gg or indeed fno- We expect the FoG effect to be in¬ 
creasingly important at smaller scales, and so these mea¬ 
surements may be at risk of misestimation. We therefore 
impose a conservative cut on the measurements, exclud¬ 
ing <7 cu t < 40 h -1 Mpc and s cut < 50 h -1 Mpc [12], In¬ 
deed, [12] showed that cosmological parameter bias be¬ 
gan to occur at smaller scales. This reduces the number 
of measurement bins in a and n to Wins = 163. 


B. Tests of theoretical templates 

When the conservative cut-off scales of (j C ut = 
40 h -1 Mpc and s cu t = 50 h~ x Mpc are used for the anal¬ 
ysis, the effective range of scale in Fourier space becomes 
k < 0.1 Mpc -1 . The power spectra of ACDM and f(R) 
gravity models are presented in this range of scale in 
Fig. 1. For the clustering scales considered in this likeli¬ 
hood analysis, there are no deviations from ACDM. This 
implies that f(R) gravity models with log \fm\ < — 6 are 
effectively equivalent to ACDM in this analysis. We take 
a uniform prior on log \Jro\ between —7 and —3. 

We first test our pipeline of analysis by checking 
whether it is possible to recover the ACDM limit 
l°g|/iJo| ^ “6 using the mock catalogues based on 
ACDM. We use the 611 CMASS mock catalogues to 
measure central values of £ (a, 7r) and fit our theoretical 
f(R ) templates to the observed correlation function. The 
measured likelihood function of log \fm\ is presented as 
a blue dotted curve in the right panel of Fig. 4. The 
best fit log | Wo I indeed lies within the ACDM limit of 
l°g l/.Ro| ^ — 6. There are no mock galaxy catalogues 
based on f(R) gravity available so we are not able to 
fully test our theoretical templates away from the ACDM 
limit. The perturbation theory predictions for the red- 
shift space power spectrum in Fourier space were, how- 














ever, tested against N-body simulations for f(R) gravity 
models and it was shown that the perturbation theory 
based template was able to reproduce the input value of 
log | fno | of the simulation in an unbiased way for k < O.f 
Mpc -1 [38], 


C. Constraints on f(R) gravity 


We now present the results for constraints on \fm\, 
summarized in Table I. The Markov Chain Monte Carlo 
(MCMC) method is used to sample a probability distri¬ 
bution. The normalised distribution function of C/C max 
of the chain is given in terms of log \/ro\ in the top panel 
of Fig. 4, and the Ay 2 is shown in the bottom panel. 

We begin by analysing our results while fixing Gq 
to the Planck ACDM model prediction of Gq = 0.46. 
The measured log|/# 0 | is then log|/# 0 | < —3.87 and 
l°g I/bo I < —3.6 at 68% and 95% confidence upper lim¬ 
its respectively. The lower index denotes the prior is 
hit - recall the ACDM limit of log |/jjoI — 6, and the 
upper indices denote 68% and 95% confidence levels of 
log | fa o | represented by unbracketed and bracketed num¬ 
bers respectively. The distribution function and Ay 2 are 
presented by blue dashed curve in the left panel of Fig. 4. 

The measured \ogfm after marginalising Gq and all 
other parameters is —4.5^^ ( ' 1 ' 4 ' ) as presented in Table I. 
The best fit \fm\ is 3.2 x 10 -5 , and the upper bounds are 
|/fio| = 3.0 x 10 -4 and 8.0 x 10 -4 at 68% and 95% con¬ 
fidence levels, respectively. The observed £(<r, 7 r) is pre¬ 
sented as blue filled contours in Fig 3, and £(<r, n) of the 
best fit f(R) gravity model is represented as black dashed 
contours, which show a slightly better fit than £(<r, n) of 
ACDM represented as black solid contours. The black 
solid curves in the left panel of Fig. 4 represents the dis¬ 
tribution function (upper) and Ay 2 (lower). Note that 
log |//jo | < —6, i.e. ACDM, is within Ay 2 < 1 of our 
best fit value and so is wholly allowed by the current 
data. When |/ fl0 | is larger than 10 -4 , £(ct, tt) becomes 
significantly different from the ACDM prediction. The 
black dotted contours in Fig. 3 represent £(<r, 7r) with 
|/fio| = 3.0 x 10 -4 at 68% confidence bound. The differ¬ 
ence from the measured £(er, n) can be observed. 

In order to check the validity of our constraints, we test 
the effect of a different initial broadband shape of power 
spectra by using the WMAP9 best fit model as fiducial 
instead of Planck 2013. Although there is no significant 
difference within statistical allowances, the predicted ini¬ 
tial broadband shape is slightly different in both cases. 
We would like to test whether this can cause a shift in 
the an apparent fno- The black dashed curve in the 
right panel of Fig. 4 represents the likelihood function 
and Ay 2 . The measured log |/jjo| using WMAP9 prior is 
—4.6^* r 1 j ( ' 1 ' 2 \ which is consistent with the result using 
Planck 2013 prior. Thus our results do not appear to be 
biased by the CMB priors. 

Finally, since we measure Jrq through its scale depen¬ 


dent effects, we consider a non-trivial scale dependent 
bias on our measurement. We model the bias as 


b(k) = b 0 


1 + A 2 k 2 
1 + Aik 


(23) 


This simple bias model, when implemented in the for¬ 
mula in Eq. (17), can nicely fit the N-body data of 
the halo power spectrum and extract correctly the lin¬ 
ear growth rate on scales very relevant to the current 
analysis [56, 57]. We find that the scale dependent 
bias parameters are determined as A 2 = —0.8±4 and 
A\ = —0.2±1.2. Given the range of errors on these pa¬ 
rameters and the range of wavenumbers k < 0.1 Mpc -1 , 
the scale dependence of the bias is poorly determined. 
Since our established cutoff on small scales removes where 
scale dependent bias is expected to be most significant, 
all other measurements, including /ro, remain similar to 
the results with scale independent bias. Thus the mea¬ 
surements shown in Table I are not modified much by a 
possible existence of this scale dependent galaxy bias. 


D. Tests of ACDM model 

Considering the constraints on the other fitting pa¬ 
rameters, we can see if there is any evidence of deviations 
from ACDM model. Information on the late-time cosmo¬ 
logical expansion is encoded in the distances along and 
transverse to the line of sight. The distance measure¬ 
ments are not affected by the scale dependent growth 
functions in f(R) gravity models. In section III, we var¬ 
ied the scale dependent growth function D / and the 
growth rate Dq and found that the BAO ring is in¬ 
variant (see Fig. 2). As a consequence, there is little 
change in the measured Da and H~ l from our previous 
result where scale independent growth functions are as¬ 
sumed [13]. The measured Da is Da = 1428.0/28.1 Mpc, 
and the measured H~ 1 is H~ x = 3157.5/[®g;® Mpc. The 
prediction of the ACDM model is indicated by the x in 
Fig. 5, and the measured values, with uncertainties, are 
presented as filled contours in the top-left panel. The 
ACDM prediction is within the 1 -a confidence level. The 
measured distances are consistent with the results using 
ACDM templates presented in the fifth column of Table I. 

Although we vary galaxy bias b in the fitting proce¬ 
dure, it is degenerate with the variation of Gg on quasi- 
linear scales. Thus as previously stated we can regard 
the combined measurement of b and Gg as a probe of Gg , 
whose measured value is 1.13/ 0 ;q 8 . Again there is little 
difference from the previous result of Gb using ACDM 
templates presented in the fifth column of Table I. 

When we fit the measured £(er, 7r), the amplitudes of 
the growth function and growth rate are determined at 
an effective scale k* ~ 0.07 h Mpc -1 in our analysis. In 
f(R) gravity the growth functions are enhanced in a scale 
dependent manner. Therefore if the amplitude of the 
growth function and growth rate at k* are tuned to be 
nearly the same as the ACDM model (e.g. to agree with 
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Parameters 

Planck 

Measurements 

Fixed G© = 0.46 

Fixed \f R0 \ = 0 

D a (Mpc) 

1397.5 

I — 1 

to 

0 ° 

o 

1 + 

to w 

00 o 

1424.3±i^ 

1422.6± 2 !;? 

H~ 1 (Mpc) 

3240.1 

3157.5^76;® 

3179 3+ 198 - 1 

ol i »-o _ 185 6 

OQi o 0 + 200.3 
OZiO.O _ 173 3 

G b 

- 

1 1 q+0-09 
x '- lo -0.08 

i io + 0 - 06 
1 - lu -0.06 

i i c+0.09 
i - lo -0.08 

G& 

0.46 

0.39l°;°9 

0.46 (fixed) 

n 41 + 0.09 
u -^ ti -0.09 

a p (Mpc) 

- 

10 - 6 ±t;| 

ll-7±l;? 

Q o+5.3 

k>g|/flo| 

- 

—4.5+ 0 : 96 / 1 - 4 ) 

prior 

_ 4 5 +0.63(0.90) 
[priori 

- 

Xmin 

- 

126.9 

127.6 

127.7 


TABLE I: The Planck best-fit ACDM predictions and the measured values with their l-cr confidence level errors are shown for 
the parameters Da, H~ x , Gb, G©, <J V , log \fm\- The Planck ACDM data does not predict the phenomenological parameters Gb 
and <T P . The last two columns show the measured values in the test cases when G© is fixed to the Planck prediction 0.46, or 
fm = 0, respectively. 


data), then the measured Gb and Gq must be smaller 
to offset the enhancement. So for appreciable log 1,/Ao I 
we expect an anti-correlation between log \fno\ and Gq. 
This is visible in the left panel of Fig. 6 for the larger 
values of log \fm\- Since the best fit log \fm\ is slightly 
larger than the LCDM bound, this may contribute to why 
the measured Gq = 0.39lQ Qg is smaller than the ACDM 
result. However, the prediction of ACDM presented by 
the x in Fig. 5 is still within the 68% CL measured Gq 
contours. 

Finally, we show the measured FoG effect in terms of 
(Tp in the right panel of Fig. 6. Again we do not find any 
significant change from the previous result using ACDM 
templates. 


IV. CONCLUSION 

We have investigated the effect of f(R) gravity on the 
redshift space correlation function £(<r, ir) using a consis¬ 
tent modified gravity perturbation theory template that 
also includes some nonlinear and chameleon screening ef¬ 
fects. The scale dependent growth functions and Dq 
alter at late times the initial broadband shape probed 
by the CMB experiments, but the comoving scale of the 
sound horizon is not affected, and we find the measured 
BAO ring remains the same as in ACDM and so the 
distances Da and iL _1 do also. The effect of modified 
gravity would be most evident in the redshift space dis¬ 
tortions and hence anisotropy of the clustering. 

We demonstrated how £(cr, 7r) alters with increments 
of D'g and Dq for small and large /rq values in Fig. 2. 
This shifts the correlation amplitude at the BAO scale 
and other scales compared with the scale independent 
growth functions case. At small \/ro\ = 3.2 x 10 -5 , 
there is no significant anisotropic variation with ADq; 
the peak locations move nearly coherently anti-clockwise 
along the BAO ring. But when \fno\ increases and it be¬ 
comes larger than \fm\ = 10” 4 , the effect of anisotropic 
amplification dominates and £(o\ 7r) departs from ACDM. 


Using the BOSS DR11 data we measure the anisotropic 
£(cr, 7r), presented as blue filled contours in Fig. 3, and 
find a best fit \Jro\ = 3.2 x 10 -5 , with an upper bound 
\fm\ < 3.0 x 10 -4 at the 68% confidence level; the lower 
68 % confidence level is consistent with ACDM. We also 
note an anticorrelation between /rq and the scale inde¬ 
pendent growth rate Gq that could reduce the measured 
value of Gq despite the enhancement of growth due to 
f{R) gravity. 

We tested our data analysis pipeline by applying it to 
ACDM mock catalogues, and we reproduced the ACDM 
bound of \Jro\ < 10~ 6 . We also tested the robustness of 
our results against the change of the CMB prior. There 
is a slight difference in the measured early broadband 
shape of spectra between WMAP9 and Planck 2013, but 
the measured \Jro\ is independent of this difference. In 
addition, we tested the effect of the scale dependent bias 
on the measured | /rq \. Again our results were shown to 
be insensitive to a possible scale dependent galaxy bias. 

Although the best fit value of |/ho| is away from ACDM 
bound, the deviation is insignificant with Ay 2 < 1, and 
it is indistinguishable from ACDM model. The compar¬ 
ison of the measured minimum y 2 nin when marginalising 
or fixing logl/Aol to zero, and marginalising or fixing 
the scale independent growth rate to the Planck best 
fit, are presented in last row of Table I. With future 
wide, deep spectroscopy experiments such as Dark En¬ 
ergy Spectroscopy Instrument (DESI) we will be able 
to use the redshift space correlation function £ (er, 7r) to 
probe f(R) gravity below \/ro\ < 10 -4 - and other scale 
dependent modified gravity - and test general relativity 
more stringently. 
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Appendix A: Power spectrum and bispectrum 
calculations in RegPT 


ized treatment is constructed with the help of the stan¬ 
dard PT kernels [59], showing that the proposed scheme 
can be used to give a percent-level prediction of the power 
spectrum and correlation function in the weakly nonlin¬ 
ear regime in both real and redshift spaces [42, 44]. 

Let us define a two-component multiplet, 4' a (fc;t) = 
(6(k\t), 0(fc; t)). Then, the power spectra of \I/ a valid at 
one-loop order are expressed as 


P ab (k-t)=T^(k-,t)T^\k-t)P 0 (k) 



x P 0 (q)P 0 (\k - q\), (Al) 

where Pq is the power spectrum of initial density field 
So- The functions F„ are the multipoint propagators. 
For the one-loop calculation, the relevant expressions for 
the regularized propagators, which reproduce both the 
resummed behavior at high-/c and standard PT results 
at low-fc, are given by 


In this appendix, based on the RegPT treatment, 
we summarize the expressions for the basic ingredients 
needed to compute the redshift-space power spectrum 
in f(R) gravity model [Eq. (17)]. The RegPT scheme 
is based on a multipoint propagator expansion, and with 
this scheme any statistical quantities consisting of density 
and velocity fields are built up with multipoint propaga¬ 
tors, in which nonperturbative properties of gravitational 
growth are wholly encapsulated [58]. Making use of the 
analytic properties of the propagators, a novel regular- 


rW(M)= k ( 1 ) (M) i + 


+3 


d 3 q 

( 2 ^ 


Fa 3 \ k ,q,-q',t)Po{q) 


,-k 2 P d /2, 


\A2) 


Fi 2) (<7, k — q-t) = F&(q, k-q;t) e~ k ^ 2 , (A3) 


(n) 

where the functions FJ are the standard PT kernels, 
sometimes written as Fa™' 1 = (F n , G n ) [43]. Note that 
the leading-order kernel F^ is related to the linear 
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growth factor through Fa ' = (-D+, -D?). The quantity 
CTd is the dispersion of the linear displacement field given 
by 

°d = J ^ p o( q ) { D +( q ',t )} 2 ■ ( A4 ) 

The above expressions are used to compute the power 
spectra Pgg, Pse> , and Poe that explicitly appear in 
Eq. (17). We also need to evaluate the A and B terms, 
which implicitly depend on the power spectrum and bis- 
pectrum [see Eqs. (19) and (20)]. Since these terms are 
regarded as next-to-leading order, the tree-level calcula¬ 
tion is sufficient for a consistent one-loop calculation of 
Eq. (17). Thus, the power spectrum and bispectrum in 
the A and B terms are evaluated with 

Pab, tree(fc) = T« (k)T^ (k) P 0 (k), (A5) 

B a bc,tree(ki,k 2 ,k 3 ) = 2T^\k 2 ,k 3 )T^(k 2 )T^(k 3 ) 

x Po(k 2 )Po(k 3 ) + (cyc.perm) (A6) 

The propagators T^ and T ^ in the above are also eval¬ 
uated with the tree-level expressions: 

rSL(*i.-" ,fc„;f)=pW(fci,--- ,fc„;t)e- fc2 ^/ 2 , 

(A7) 


where k = |fei... n |. 


Note finally that while the power spectrum and bis- 
pectrum given above are quite general and valid in any 
model of modified gravity, the expressions for the prop¬ 
agators in RegPT might receive some corrections. This 
issue has been discussed in detail in Ref. [10]. However, 
it has been found that within the f(R ) gravity model, 
any corrections due to the modification of gravity are 
very small and can be neglected at the large scales of 
our interest. Hence, all the expressions given above are 
basically the same as those in GR, except for the stan¬ 
dard PT kernels, ■ Unlike the GR case, the time and 
scale dependence of Fa are not separately treated in 
f{R) gravity model. The functional form of Fa' l> is thus 
non-trivial, even if the deviation of gravity from GR is 
small. In this paper, on the basis of Eqs. (2), (3), (6) 
and (7), we derive the equations that govern the stan¬ 
dard PT kernels. Solving these equations numerically, 
we construct the kernels that are tabulated as function 
of wave vectors [60]. 
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